Phase transitions in a mechanical system coupled to 
Glauber spins 



A Prados\ L L Bonilla^ and A Carpio^ 

^Fi'sica Teorica, Universidad de Sevilla, Apartado de Correos 1065, E-41080, Sevilla, 
Spain 

^ G. Millan Institute for Fluid Dynamics, Nanoscience and Industrial Mathematics, 
Universidad Carlos III de Madrid, 28911 Leganes, Spain 

Departamento de Matematica Aplicada, Universidad Complutense de Madrid, 
28040 Madrid, Spain 

E-mail: pradosSus . es,bonilla@ing.uc3m. es,carpio@mat .ucm. es 

Abstract. A harmonic oscillator linearly coupled with a linear chain of Ising spins 
is investigated. The N spins in the chain interact with their nearest neighbours 
with a coupling constant proportional to the oscillator position and to N~^/^, are 
in contact with a thermal bath at temperature T, and evolve under Glauber dynamics. 
The oscillator position is a stochastic process due to the oscillator-spin interaction 
which produces drastic changes in the equilibrium behaviour and the dynamics of the 
oscillator. Firstly, there is a second order phase transition at a critical temperature 
whose order parameter is the oscillator stable rest position: this position is zero above 
Tc and different from zero below Tc. This transition appears because the oscillator 
moves in an effective potential equal to the harmonic term plus the free energy of 
the spin system at fixed oscillator position. Secondly, assuming fast spin relaxation 
(compared to the oscillator natural period), the oscillator dynamical behaviour is 
described by an effective equation containing a nonlinear friction term that drives the 
oscillator towards the stable equilibrium state of the effective potential. The analytical 
results are compared with numerical simulation throughout the paper. 
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1. Introduction 

Many physical processes are interpreted in terms of an oscillator coupled to a 
thermal bath or to spin systems. Examples abound, mass spectrometry through a 
nanoelectromechanical oscillator whose resonant frequency decreases as single molecules 
are added to it [1] , a spin representing a two- level system is coupled to a boson bath (the 
spin-boson system) to analyze loss of quantum coherence due to the bath [2], a classical 
oscillator coupled to a spin causes wave function collapse thereof [3], the classical version 
of the spin-phonon system describes the collective Jahn- Teller effect |H |5], large spin 
systems (single molecule magnets or nuclear spins) are coupled to a boson bath [6] , etc. 

In this work, we consider a mechanical degree of freedom represented by a classical 
harmonic oscillator coupled to a linear chain of N Ising spins ai {i = 1, . . . , N, ai = ±1) 
in contact with a thermal bath at temperature T. The energy of the combined system 
is equal to the energy of the oscillator alone plus a coupling term proportional to the 
oscillator position and to N~^^^ J2iLi <^i<^i+i- The spins flip stochastically according to 
Glauber dynamics [7j. As a consequence of the coupling, the oscillator equations of 
motion become stochastic, and both the position and the momentum of the oscillator 
become stochastic processes. The aim of this work is to understand how the equilibrium 
and the dynamics of the oscillator is affected by the interaction with the spin system 
(and vice versa). 

The plan of the paper is as follows. The oscillator-spin model is described in 
section O Its time evolution is governed by Newton's second law for the oscillator and 
the above mentioned Glauber dynamics for the spins. They can be put together in an 
evolution equation for the joint probability density of finding at time t the oscillator at 
given values of its position and momentum and the Ising system at a given configuration. 
The canonical distribution at temperature T is the equilibrium joint probability density. 
By summing over all possible spin configurations, we obtain the equilibrium probability 
density for the oscillator. The latter is a canonical distribution with an effective potential 
energy which is the sum of the harmonic potential and the free energy of the spin chain 
for fixed oscillator position. 

Section |3] is devoted to analyzing the equilibrium configuration. By finding the 
minima of the effective potential, we show that there is a second order phase transition 
at a critical temperature Tc, with the stable rest position of the oscillator (equilibrium) as 
its order parameter. For T > Tc, the oscillator equilibrium position is the same as that of 
the uncoupled oscillator. For T < T^, two symmetric nonzero equilibrium positions issue 
forth from zero as in the diagram of a pitchfork bifurcation. These nonzero equilibrium 
positions behave as \/iV as the number of spins goes to infinity, so that the harmonic 
contribution to the energy be extensive in the thermodynamic limit. On the other hand, 
the fluctuations scale as N~^^'^ far from the critical temperature. Very close to Tc, there 
is a crossover and equilibrium fluctuations scale as A^~^/^. Although Ising spins in the 
chain are coupled to their nearest neighbours, their coupling constant is proportional 
to the oscillator position which makes their interaction effectively long range. Similar 
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hidden long-range effective correlations that enable possible Id phase transitions are 
present in biophysical systems. An example is DNA melting |8] which has been modeled 
by means of modified Ising systems [9l ^U\, different from the one considered here. 

The dynamics of the system is studied in section HI In the limit of fast relaxation 
of the spins compared to the natural period of the oscillator, there is a clear separation 
of time scales, a fast one associated to the relaxation of the spins and a slow one 
associated to the oscillator. In this regime, we find a reduced dynamics of the oscillator 
with nonlinear friction and a nonlinear force term. This nonlinear evolution equation 
is one of the main results of our paper. Basically, the spins approach their equilibrium 
distribution corresponding to the instantaneous value of the oscillator position. This 
produces the effective potential (already found in the equilibrium analysis) for the 
oscillator and gives rise to the nonlinear force term in its nonlinear evolution equation. 
On the other hand, the nonlinear friction is a purely dynamical effect that cannot be 
obtained from analyzing the equilibrium distribution of the system. This friction arises 
from the slow evolution of the oscillator resulting in a slight deviation of the Ising spins 
from its equilibrium with a fixed position of the oscillator. The friction term drives the 
system to equilibrium in the long time limit. The stationary solutions of the reduced 
dynamics coincide with the oscillator equilibrium positions at any given temperature. 
We also discuss the expected range of validity of the nonlinear dynamical equation. For 
T > Tc (section HIT]), we can linearize the oscillator reduced evolution equation about its 
stable rest state. The solutions are underdamped oscillations whose frequency decreases 
as T decreases: the oscillator is slowed down by the spins. There is a narrow region of 
overdamped oscillations for temperatures very close to T^. A similar analysis is carried 
out in section 14.21 but for T < Tc. There is also a very narrow region of overdamped 
oscillations near Tc. For lower temperatures T < T^, our theory predicts underdamped 
oscillations around one of the two nonvanishing stable equilibrium points. 

In Section [5l we compare numerical simulations for the model with the theoretical 
results and test the range of validity of the theory. The numerical simulations show 
excellent agreement with the theory for sufficiently high temperature T > Tc and for 
1000 or more spins. As T decreases towards Tc, the numerical solutions of our theory 
and the simulations show the same qualitative trends, i.e., underdamped oscillations, 
but these oscillations are shifted by some constant value. This is an effect due to 
the initial conditions as reducing the size thereof brings again quantitative agreement 
between theory and simulations. Below the critical temperature, but not very far from 
it, more spins are necessary to attain good average values and our theory still gives an 
adequate description of the dynamical evolution of the system. As the temperature is 
further lowered, we again need fewer spins to attain good averages over spin indices 
and trajectories but there are qualitative differences between theory and simulations. 
Breakdown of the theoretical predictions is expected for sufficiently low temperatures, 
because of the divergence of the relaxation time of the spins [HI [121 [13] • Lastly, section 
E] contains final remarks and comments. 



Phase transitions in a mechanical system coupled to Glauber spins 



4 



2. The model 

We consider a system comprising one dimensional harmonic oscillator (mass m, 
frequency Uq, position x and momentum p) and ^ 1 internal degrees of freedom 
modeled by Ising spins (cxj = ±1, i = 1,...,N) in contact with a heat bath at 
temperature T. The system has an energy 

'H{x,p,a) = 'Ho{x,p) +'Hintix,cr) (la) 

N 

'Hint(a;,cr) = - /ix^o-iCJi+i , (Ic) 

i=l 

in which 'Ho{x,p) and Tii^tix, cr) are the energy of the uncoupled oscillator and the 
interaction energy between the oscillator and the spins, respectively. The latter can 
also be understood as a nearest neighbour interaction between the spins with a coupling 
constant Jes which is proportional to the oscillator position x, 

JeS = fix . (2) 

The parameter /i measures the strength of the coupling between the oscillator and the 
Ising system. Because of the sum over spins in ( |Tc|) . /i should decrease with for the 
system to have a well defined behaviour in the limit N ^ oo. We will show later that 
/X = /io/ as mentioned in section [H Alternatively, the Hamiltonian ( II ap can also be 
written as 

n{x,p,(T) = — + V{x,(t) , (3a) 
1 ^ 

V{x,cr) = -mulx^ - /xx (JifTj+i , (36) 
^ i=i 

where V(a;, cr) is the total potential acting on the oscillator. 

The dynamics of the system is governed by Hamilton's equations of motion for the 
oscillator, 

X = — (3a) 
m 

P = - t;^ = -mujlx + yu ^ o-.o-i+i , (36) 

i=l 



or, equivalently, 

N 

m 



X + ojIx = — a-iCXi+i , (3 c) 

i=l 

and by an appropriate stochastic dynamics for the spins (because they are in contact 
with a heat bath at temperature T). For the sake of simplicity, the spins will be 
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assumed to evolve with Glauber-like one spin flip dynamics. At any time t, the system 
may experience a transition from {x,p, a) to (x, p, Ricr) with a rate given by [7] 



Wi{cr\x,p) 



a 



1 ~ + CTj+ly 



where RiCr is the configuration obtained from cr by rotating the i-ih. spin. Here 



7 = tanh 



2 J, 



cff 



tanh 



(4) 



(5) 



A;^ is the Boltzmann constant and T is the temperature of the system. The quantity a 
determines the characteristic attempt rate for the transitions in the Ising system. 

In this way, the joint probability V{x,p, cr, t) of finding the oscillator with position 
X and momentum p, and the spins in a configuration cr = {ai, (J2, . . . , Cn} at time t 
obeys the Liouville-master equation 



dtV{x, p, cr, t) + ^dj^{x, p, cr,t) + | -mcUgX + /i ^ o-jCTj+i | dpV{x, p, cr, t) 

N 



i=l 



J2 [Wi{Ri(T\x, p)V{x, p, RiO-, t) - Wi{(T\x, p)V{x, p, a, t)] 



i=l 



The equilibrium solution of this equation is the canonical distribution 



where Z is the partition function 

/ + 00 /■ + 00 

dx / dp < 
-OO J —CO _ 



-/3W(x,p,o-) 



(6) 



(7) 



(8) 



and (3 = {ksT) ^. Since we are mainly interested in the behaviour of the oscillator, it 
will be useful to consider the marginal probability Vcci{x,p) 

1 
Z 



(9) 



where 



ZisUx) = e-^-^i--g 



2 cosh 



cff 



1 N 



(10) 



is the partition function of a Id nearest neighbour Ising model with coupling constant 
Jeff, which depends on x as given by ([2]), and J-ising(a;) the corresponding free energy. 
Therefore, Veq{x,p) is readily rewritten. 



'Peqix,p) = — exp <j -13 



2m 



with 



'^esix) = -mUQX + Jlsing(x 



-muolx'^ — NksT 



In cosh 



[IX 



+ ln2 



(11) 

(12) 
(13) 
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Equation ( ITT]) suggests that Ves is the effective potential acting on the oscillator due to 
its coupling to the Glauber spins. This point will be confirmed when the dynamics 
be analyzed in section HI 



2.1. Orders of magnitude and nondimensional equations 



X 


P 


V 




t 


e 


e 


liN 






a 


j_ 




T muigkBT 






11^ 




a 





Table 1. Nondimensional units and parameters. 



It is convenient to render our equations dimensionless before we proceed with their 
analysis. To do this, we can start with Eq. ( !3cj) . The two terms in its left hand side 
have the same order if we adopt t* = u^t as a nondimensional time. The spins (Tj are 
either +1 or -1, and therefore its right hand side (the forcing term) is, at most, fiN/m. 
Adopting this value as an order of magnitude of the forcing term, it is of the same order 
of magnitude as any of the terms in the left side of fl3cj) provided x has an order of 
magnitude [x] = fiN/{mu)Q). The normalization condition 



oo 



J2 dx / dpV{x,p,a,t) =1, (14) 



yields 



[x] [p] muolx]'^ p?N'^ 
Lastly, the argument of the coefficient in Eq. ([5]) has order of magnitude 

^ [x] fi^N 



A;bT mujQkBT T' 



where 



Tc = (15) 

is a critical temperature whose role we will unveil later in the paper. 

Thus we can define nondimensional variables according to x* = x/[x], t* = t/[t], 
. . . , where the units [x], [t], ... are as defined in Table [H Inserting these nondimensional 
variables in Equations ( l3cj) . (jl]), ([5]) and (Q, and dropping the asterisks in the result (so 
as not to clutter our formulas), we obtain the following nondimensional equations 



=1 



N 

J2 [W,{R,(T\x,p)V{x,p,R,(T,t) - Wi{(T\x,p)V{x,p,(T,t)] 

i=l 
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V{x,p,(T,t), 



dt+pd^+ ^ (r^(r^+l - dp 



1 ^fx) 

Wi{cr\x,p) = —o-i{ai^i + o-,;+i), 



7(x) = tanh ( — J , e 



2x\ loq 

In nondimensional units, the equilibrium distributions ([7]) and f lTT]) are 



Peq(a;,p, cr) = — exp 



AT 

T 



n{x,p, cr) 



p + X X 

— ^ ]V^''''^'+'' 

i=l 



and 



Peq(a;,p) = - exp 



9 V 2 



a; 



+ Vcff(x) 



In cosh ( — ) + In 2 
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(17) 

(18) 
(19) 

(20) 
(21) 

(22) 
(23) 



respectively. 



3. Equilibrium points and phase transition 

The maxima of Veq{x,p) determine the most likely position and momentum of the 
oscillator coupled to the Ising system, (xeq,Peq), when the total system is at equilibrium. 
These most likely values will be called macroscopic equilibrium values following van 
Kampen's terminology [Hj. As N oo, the equilibrium mean values of x and p 



coincide with Xeq and Peq, respectively, whereas the corresponding variances tend to 
zero. Similarly, the equilibrium average value of an smooth function f{x,p) tends to its 
macroscopic value: {f{x,p)) ~ f(xeq,Peq) as — )■ oo. Thus a macroscopic quantity has 
negligible fluctuations in the limit of infinitely many oscillators. Let us now calculate 
Xeq and and the corresponding variances. First, Pcq = and the oscillator is at rest 
in equilibrium, as expected. Second, the oscillator macroscopic equilibrium positions 
are given by the solutions of the equation 
dVeff(a;) 



dx 



0. 



I.e., 



X 



eq 



tanh 



eq 



9 



0. 



(24) 



(25) 



Clearly Xeq = is always a solution for any value of 9. It is the only solution for 9 > 1, 
it corresponds to a maximum of Veq and is therefore stable. At 9 = 1 two new stable 
equilibria issue from x^q = and exist for 9 < 1. Note that in dimensional units, 
9 = 9c = I corresponds to T = Tc, the critical temperature defined in ( fTSl) . Besides, Tc 
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should be independent of in the large N limit. This gives the scaling of /i with 
mentioned in the Introduction, 



where Hq is independent of N. Therefore, 



making use of f|T5|) . 

As — )■ 1~, we find 



(26) 



(27) 



i.e. the usual scaling at pitchfork bifurcations, 
fl23l) is continuous at 6* = 1, 



v: 



:^ + ^ln2~- 



1 



39 (\ 



(28) 

The effective potential 



eq 



(29) 



as 6* —7- 1^. Then the derivative of V^g with respect to T is also continuous at = 1. We 
have found a second order, continuous, phase transition with classical critical exponents. 
The equilibrium position x can then be considered the order parameter of the transition: 
its macroscopic value vanishes for > 1 and is non-zero for < 1. 

Why does this second order transition appear? At first, it seems surprising to find 
it in a Id model with short-ranged interactions. In order to understand the physical 
reason for this behaviour, let us calculate the equilibrium probability Peq(c) of finding 
the spins in configuration o", regardless of the values of x and p. We shall integrate the 
probability density ( |20l) . written as 

X-^{(T)) — 



Veq{x,p,(7 

where 



^exp 



A^ 

T 



p^ 1 
— + - 

2 2 



over X and p, with the result 

1 



exp 



exp 



29N 



(30) 



(31) 



(32) 



Here Si = ajCTj+i are new effective spin variables and Z^f is the appropriate normalization 
constant. Then ( 132|) corresponds to the equilibrium probability of a mean field Ising 
model. Each spin Sj is coupled to the global mean field ip = Sj/N and an effective 
long range interaction appears in the model. It is a well-known result that the Id 
mean field Ising model has a second order phase transition at a finite temperature |15j . 
The macroscopic, most probable, value of ip is given by solutions of the trascendental 
equation [151 fT6] 



tanh 



9 



(33) 
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There appears a second order transition at a critical temperature 6 = 1, which is the 
same one appearing in Eq. fl28|) . The origin of this transition is the hidden long-range 
effective coupling between spins ( l32l) which is produced by the coupling of the Glauber 
spins to the oscillator. Similar hidden long-range effective correlations that enable 
possible Id phase transitions are present in biophysical systems. An example is DNA 
melting |8] which has been modeled by means of modified Ising systems |9l[Tn], different 
from the one considered here. For DNA melting, the order of the phase transition has 
not yet been well estabhshed: depending on models and conditions, it has been predicted 
to be first order [T7J, second order pTS], or even higher [19]. 

The fluctuations of the order parameter x can be analyzed from the equilibrium 
distribution fl22l) in the limit — ?■ oo. The average of any function of x can be calculated 
by using the Laplace method in integrals involving ( l22l) . which leads to expanding the 
effective potential ( l23l) around the macroscopic value Xgq- The result is 



where 



N 



1 



N . 
—bj- 
2 



X 2^eq ) ) 



1-^(1-., 



eqj 



(34) 



(35) 



w is a new dimensionless frequency. Therefore, for u ^ {6 ^ 1), the fluctuations of 
X are Gaussian because higher order terms vanish as A^ — t- oo. The average value of x 
equals Xeq, as expected, and its variance is 

'1 



<7^ 



((x Xgq) ) 



9 



eq 



X 



eqy 



(36) 



A^ 



(Vefr)eq " V, 



(37) 



Nu^ N 

which vanishes as A^~^. Similarly, 0"^ = (p^)cq = 0/N. Since Xgq is of order one for A^ ^ 1 
and 6^1, the fluctuation of x is much smaller than its average value, which is the 
expected behaviour of a macroscopic variable. The average value of Ves at equilibrium 
verifies 

— U ((X-Xeq) ) = - 

for 9^1. The term coming from Gaussian fluctuations is subdominant in the 
thermodynamic limit as compared to the extensive macroscopic contribution A^V^^. 

On the other hand, w — )■ and therefore ax in ( 136|) diverges as 6* — )■ 1: fluctuation 
divergence is connected to the vanishing of the renormalized frequency uj. A very 
large value of A^, diverging for 6* — )■ 1, has to be considered in order to be in the 
"thermodynamic limit" for the oscillator position, where x is approximately equal to 
its most probable value Xeq and its fluctuactions can be neglected. What happens for 
9 = 11 The first three differentials of the effective potential vanish at 6* = 1 whereas 
/dx^ = 2/9^. Then, as ^ 1+ 



eq 



1 ^Xj 3y I 



cq) 



+ 



N{x 



X, 



eqj 



12 



(3^ 
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and a similar expression holds as 6 ^ 1^ (replacing {0 — 1) in (155]) by 2(1 — 6')). Therefore 
the fluctuations scale is (x — x^q) oc N^^^^ as — oo if — 1| <^ N^^^'^ ^ 1 (non- 
Gaussian behaviour, the quadratic term can be neglected in comparison to the quartic 
term) and N ^^'^ a N ^/'^ <^ \9 — 1\ (Gaussian behaviour, the quartic 

term is negligible). 



4. Dynamics 



In this section we shall analyze the dynamical equations of motion. Equation (fT6ll is 
a stochastic differential equation for x because the configuration of the spin system cr 
is a stochastic process. Let us denote Cj_„ = crjCrj+„, with i = 1, . . . , N and n > 0. Of 
course, Cifi = af = 1 for all i. By averaging ( |T6l) over the joint probability V{x,p, cr,t) 
solution of the Liouville-master equation ([6]), we obtain 

+ {-) = jjY.iC..) ■ (39) 

1=1 

From the Liouville-master equation, we can derive the following system of equations for 
the spin correlations 

- 2(Ci,n) + ^(7(3;) (Ci,„-1 + + Ci-l,n+l + Q+l.n-l)) = ^ ^ , (40) 

for n > 1 and i = 1, . . . , N. Here 7(0;) and e are given by f|T9|) . The system of equations 
( 140|) must be solved with the boundary condition Ci^o = 1 and given initial conditions 
{(a,n)(t = 0);n> 1}. 

As explained in section [HI a quantity is called macroscopic if, compared to its mean, 
its fluctuations are negligible in the limit as — )■ 00. In this section, we will describe 
the mean-field (macroscopic) dynamics of our oscillator-spin system such that 



(F(x,p, 0-i(Ti+„))(t) 



l,xl,pj:Fix,p,...M,p,.,t) 



^ F{x{t),p{t),C^n{t)), (41) 
(a,(t)a,;+„(t)) = (a,„)(t)~Q,„(t), ^ = l,...,Ar, (42) 

in the limit as A^ — )■ 00 for any smooth function F{x,y,z). In these equations and 
for each time t, x{t), p{t) and Ci^n{t) are the values of x, p and of Cj^„ for which the 
probability density function V{x,p, cr, t) has a maximum. Due to translation invariance, 
in the limit as A^ — )■ cxd, the averages {Ci^n){i) are independent of i, provided the initial 
probability density is translation invariant and isotropic. Then Ci^n{t) is independent 
of i and we can write Cn(t) instead of Cj,„(t) (or {Ci^n){t)) in (l39l)-fH2l). Since af = 1, 
we have Cq = 1. Ignoring fluctuations according to (14T|) . equation (!39l) yields 

d^T - — 

d^ + 5^ = ^i, (43) 
and fl4Up simplifies to 

-2Cn + 7(J) (aZ'i + CUi) = e ^^Cn , n > 1 , (44) 
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corresponding to van Kampen's macroscopic approximation [13]. This approximation 
is equivalent to separating macroscopic and fluctuating contributions in x and Cj^„: 

X =x + 6x, (45) 

Q,n = Cn + 6a,n. (46) 

Inserting these expressions in (l40l) and neglecting all terms containing correlations, such 
as ((5x)^) or {6x6Ci^n), we obtain again dH]). The mean-field or macroscopic dynamical 
behaviour of the oscillator-spin system is found by solving the equations ( l43l) and ( l44l) 
with the boundary condition Cq = 1 and appropriate initial conditions. 

We now consider the limit e = a;o/a<^lofa very slow oscillator compared to the 
relaxation time of the Glauber spins. Setting e = in f|T7j) . we find the equilibrium 
solution of the master equation for each instantaneous value of x{t). For e ^ 1, there is 
an initial time window inside which the Cn reach their equilibrium values (corresponding 
to the equilibrium solution of the master equation with fixed x{t)) while the oscillator 
position and velocity are frozen at their initial values. After this initial layer, we can 
approximately solve pUj) by means of an expansion in powers of e: 



Cr.{t;e) = Y,Cn\t)e' + 0{e'). (47) 

This yields the boundary conditions C^^ = 1, Cq^^ = 0. Inserting (H71) into (H^ . we 
obtain the following system of equations 

7(J)(Cf^ + cSi) - 2Cf = 0,^ (48) 

7imci!l, + Ci!l)-2Ci'^ = ^, (49) 

and so on. The solutions of fH5]) and (jlSD with boundary conditions C^^^ = 1 and 
Cg^"* = are found in Appendix A They provide 



~ /x\ e 1 + tanh^ (f ) dx 

Ci = tanh - -— TT^^- 50 

\e J 201-tanh2(|) dt ^ ' 

Eq. ( 150|) comes from a "normal" solution of the system of equations [201 EI] , in which all 
the time dependence in Cn occurs through x, which evolves on the slower time scale t. 
The first term in fl50l) is the equilibrium value of Ci corresponding to the instantaneous 
oscillator position x, whereas the second term contains the (small) deviations from 
equilibrium, to order e ^ 1. The initial condition Ciit = 0) does not appear in 
expression (!50!) because the spins forget their initial conditions on a time scale (initial 
layer) much shorter than the natural period of the oscillator. Inserting (|50l) into ( H3l) . 
we get 
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Equation f lSip can be rewritten in terms of the nondimensional effective potential 
fl23|l and the friction coefficient 

l+tanh^(f) , , 

R(x) = 52 

^ ' l-tanh'd)' ^ ' 

as 

f --V^.(S)-^fl(J)f. (53) 

The equihbrium values of x can be obtained from fl^Tl) and the results of the previous 
section are recovered. Again stable equilibrium points correspond to the minima of Veff. 

Equation (15T1) . or equivalently (!53|l . is the main result of this section. It shows 
that the effect of the coupling of the oscillator with the bath of Ising spins is twofold. 
Firstly, the potential is "renormalized" to Ves, as a new force tanh(5;/6') is added to 
the harmonic interaction —x; secondly, a nonlinear friction term proportional to dx/dt 
appears. 

Which is the expected range of validity of the nonlinear equation (1511) ? For a given 
value of e, the range depends on the order of x and 6. Equation f lSTj) holds if the spins 
relax to equilibrium so fast that the oscillator position does not change. Using that 
Am = 2[1 — 7(x)] is the smallest eigenvalue of the coefficient matrix in ( l44l) \7], we find 
A„i ~ 2 for 6* x (high temperature limit) and A™, ~ e for 



1 /^2^ 
tanh I — 



|--7lnf7)- (54) 



9 4 V4. 

For x/6 satisfying f lM|) and larger values, our separation of time scales breaks down 
and we do not expect (151]) to hold. Due to the logarithmic dependence on e in ( 15^ . 
our asymptotic theory should already fail for moderate values of x/6. For instance, 
x/9 ^ 1.5 for e = 0.01. More will be said about this point in the numerical section. We 
now particularize our theory for temperatures above and below critical. 

4.1. The region 6 > 1 

In the high temperature region 6 > 1, the stable equilibrium point of the oscillator is 
Xeq = 0. This equilibrium is asymptotically stable due to the presence of the damping 
term, and therefore x{t) -C 1 for long enough times. Then f lSTj) can be approximated 
as: 

d^x e dx / 1 
dt^^2^dt ^ \9, 

which is the equation of the damped harmonic oscillator with square frequency 1 — 1/^ 
(which equals the renormalized frequency defined in (15^ for x^q = 0) and friction 
coefficient e/{26). The renormalized frequency tends to zero as ^ — ?• 1^. 
Defining the damping ratio 



+ — — + 1--U = 0, (55) 
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the underdamped, critically damped and overdamped oscillations correspond to C < 1, 
C = 1 and C > 1) respectively. Thus, a new dynamical "critical" temperature 6'^ appears 
for > 1, defined by the condition C = 1- According to ( l56l) . this occurs for 

Then, in the limit e = cuo/a ^ 1 we are analyzing, is very close to the critical 
temperature 1 and the region of overdamped oscillations is very narrow: its width is of 
the order of e^. 



4-2. The region 9 < 1 

In this region, the equilibrium position of the oscillator is given by the nonvanishing 
solutions of (133|) . ±Xeq- Thus the dynamics will be governed by the nonlinear equations 
If Xi(t) is one solution evolving towards Xeq as t — )■ oo, — is also a solution 
which evolves towards — Xeq- This is not in contradiction with the /mear Liouville- master 
equation having a unique equilibrium distribution fl2CT]) -f l?I]) . In fact, the multiplicity 
of (macroscopic) equilibrium solutions corresponding to extrema of the equilibrium 
distribution is a direct consequence of the nonlinearity of the macroscopic equation, 
which is compatible with the linearity of the master equation [H]. 
Near the stable equilibrium points, 

X = Xeq + ^, (58) 

with ^ ^ 1, we linearize (151]) or (!53l) . thereby obtaining 

§ + r_| + .^^^0. (59) 

Here, 

c^! = V:'g(xeq) = 1 - (60a) 
=^ R(x,.) = (606) 

The frequency U- is equal to the renormalized frequency u introduced in ([35]), 
particularized for 6 < 1. As in the case ^ > 1, we find the equations of a damped 
harmonic oscillator, but with different friction coefficient and frequency. As — ?■ 1^, we 

get 

^2 ^ 2 (1 - ^) , (60c) 

which also tends to zero. 

The analysis of the overdamped, critically damped and underdamped oscillation 
regions is completely analogous to the case 6 > 1. We find a new temperature 6^, 
for which the oscillations are critically damped. Thus for d < the oscillations are 
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underdamped, while ioi 6^ < 6 < 1 they are overdamped. The critical dynamical 
temperature is determined by 

r 

which gives, after some calculation, 

= 1 - 1^ + O(e^) . (62) 

Again, the region of overdamped oscillations below = 1 is very narrow. It should be 
noted that another region of overdamped oscillations is predicted by ( 159|) - (]606|) for very 



low temperatures, as the friction coefficient r_ formally diverges for T — )■ or x^q 1- 
Nevertheless, we will not investigate this region because it lies outside the range of 
validity of our dynamical equation ( 15T|) . as we discuss in relation with the numerical 
results in the next section. 

5. Numerical results 

In order to test our theoretical predictions, we have carried out numerical simulations 
of the stochastic process corresponding to the dimensionless Liouville-master equation 
( IT7|) -( [T8|) . Equivalently, we have to integrate numerically the oscillator equation (fT6|) and 
the Glauber evolution equations for the spins given by the transition probabilities ( ITSl) . 
Using the initial probability distribution V{x,p, (t,0), we generate initial conditions 
(a;(0; z/),p(0; z/), cr(0; z/)) for Nt trajectories u {v = 1, • • • , Nt). The oscillator position 
and momentum and the spin configuration of a given trajectory u at time t are denoted 
by x{t] u), pit; u) and ait; u), respectively. For a given trajectory at time t, we choose 
at random one spin cxj and fiip it with probability Wi{(T\x{t] u) , p{t; u)) < 1 at time 
t + At (in our dimensionless time scale with time unit At = e/N, according to 

the Metropolis algorithm for the master equation [221 123]; iii dimensional units, we have 
At = e/{Nuo) = (Na)^^ ). The oscillator position and momentum are also updated by 

x(t + -^;u^ = x{t; + P{t; ^) , (63) 



P{^ + J;^'^) =P(^!^) + ^ -x{t;u) + ^^ai{t;u)a,+i{t;u) 



(64) 



Up to a certain time to, each trajectory u is obtained by iterating this procedure 
toNe'^ times. Afterwards, the numerical averages over u (trajectories) give the averages 
with the probability distribution V{x,p, cr). In particular, we get the self-averaging 
properties, 

^ Nt 
^ Nt N 

(CiW) = ^$^$^a,(t;z/)a,+i(t;z/) ^Ci(t), (66) 

^ u=l i=l 
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as — i- oo and oo. The simulation results {x{t)) and {Ci{t)) should therefore 

approach the macroscopic values x{t) and Ci(t), respectively, for sufficiently large N and 
Nt- The numbers of spins and of trajectories Nt needed to get good approximations 
to X and Ci in (l65ll and ( !66|) are related to the amplitude of the averaged trajectories. 
As this amplitude decreases, the number of particles and trajectories must be increased. 
When 5;(0) = 0(1) and x{0) = 0(1), good averages are obtained with > 10^ and 
Nj- > 10^. We have used N = 10^ and Nt = 10^ in our numerical simulations although 
we have observed that, depending on the initial values of x and p, we can take smaller N 
and Nt without losing accuracy. Note that an order-one initial dimensionless position 
X of the oscillator corresponds to a dimensional position x of order a/ZV. This means 
that the oscillator energy is comparable to the energy of the spin system, which is also 
of order N. 

Our theory is expected to provide a good description of the numerical curves if 
e ^ 1 and x{0)/9 is not too large (high temperature). As the temperature decreases for 
fixed x(0), the characteristic relaxation time of the Ising system increases and becomes 
comparable to the oscillator period when x{0)/6 satisfies fl54|) . For lower temperature, 
we expect our theory to break down. Let us check this from the results of the numerical 
simulations. 

In figure [U^a), we show the time evolution of the oscillator for one high value of 
the temperature, namely 6 = 4, corresponding to the underdamped region. We have 
chosen initial conditions so that the spin system is initially in a completely random 
state (therefore Ci(0) = 0) and, for the oscillator (5;(0), 5;(0)) = (1,1). The numerical 
curves have been obtained with N = 10^ spins and averaged over Nt = 10^ trajectories. 
The theoretical predictions based on both the nonlinear evolution equation (ISTl) and the 
linear approximation fl55|) (which are almost indistinguishable) show excellent agreement 
with the simulations. For the plotted values, tanh(5;/6') c^x/9 thereby justifying the use 
of the linear approximation ( l55l) . Similar behaviour is obtained for different values of 
x(0) and x(0), provided x/6 1 for all times. For larger x(0), the nonlinear evolution 
equation ( I5T1) describes well the dynamics but there is a initial time window for which x/9 
is not small and the linear approximation is valid. Since given sufficient time, x — )■ for 
any initial condition, the linear equation ( l55l) always provides a good approximation of 
the dynamics for long enough times. Similarly, for the same initial conditions as in Figure 
[U^a), Figure [D^b) shows that the nonlinear equation (ISTj) gives a good approximation 
of the simulations for a lower temperature 6 = 2 but the linear equation (l55ll does not: 
tanh(x/6') ^ 5;/6' no longer holds. 

For ^ = 1.1 (a value closer to the critical temperature ^ = 1) and the same initial 
conditions. Figure [2]^a) shows that even the predictions based on the nonlinear equation 
( 15T|) fail to approximate the simulation results. This could have been foreseen because 
Eq. (jSlj) gives 6 ~ 1.44 (for an initial a; ~ 1) as the limiting temperature above which 
the nonlinear equation holds. Figure [2]^b) shows that the difference between the average 
value Cf™ in the simulation and the theoretical prediction ( l50l) . is rather larger 
than the theoretical error of order = 0.0625 for times t < 10. 




Figure 1. Averaged trajectories {x{t)) = x{t) (circles) versus nonlinear (solid blue 
line) and linear (dot-dashed green line) predictions for initial data x{0) — 1, i(0) = 1 
and (a) = 4, (b) 6* = 2. Other parameter values are e = 0.25, N = 10"^, Nt = 10^. 



Selecting again the completely random state (Cf™(0)=0) as initial condition of the 
Ising spins, we have considered smaller values of x(0) and x(0). They are such that 
initially the rhs of ([50]) is of order and Cf'^(O) - ^^(0) = 0{e^). Therefore, (jSOD is 
initially valid and the possible departure from ( ISTl) cannot be a transient effect, due to 
"inadequate" initial conditions. Equation (|50|1 suggests immediately the choice 

x(0) = e^ x(0) = e. (67) 

With this choice of initial conditions, the nonlinear equation (|50l) is still a good 
approximation for the dynamical behaviour of the oscillator at temperature 6 = 1.1, 
as shown by figure [3l Interestingly, this choice corresponds to variables of the order of 
unity for a different, alternative, nondimensionalization of the variables, 

z = x/e^, T = t/e, (68) 
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Figure 2. (a) Comparison between tlie averaged trajectories (red circles), the 
nonlinear equation ([?T|) (solid blue line) and the linear equation (1551) (dashed green 
line) for 9 = 1.1, near but above the transition temperature, (b) Difference between 
the average value Cf^™ in the simulation and the theoretical prediction Cf\ ([SO]) . 

1| . . . . 1 
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Figure 3. Averaged trajectories {x{t)) — x{t) (circles) versus nonlinear (solid blue 
line) and linear (dot-dashed green line) predictions for initial data x(0) = e^, x{0) = e 
and 9 = 1.1. Other parameter values are as in Figure [T] 



in which the spin relaxation time 1/a is selected as the unit of time instead of I/cuq as in 
table [TJ This choice is the natural one to monitor spin relaxation at high temperature. In 
Fig. m we observe that the linear approximation breaks down for temperatures closer or 
equal to the critical value 6 = 1 but the nonlinear equation is still a good approximation 
of the simulation values. In particular, this is the situation in the overdamped region 
9^ < 9 < 9^ , where 9"^ are given by ( 1571) and (1621) . respectively. For the value e = 0.25, 
9'^ = 0.998 and 6'j" = 1.004. It is a very narrow region, being its width of order e^. 
Therefore, the evolution of the oscillator is almost indistinguishable from the critical 
temperature behaviour, shown in figure HJ^b). 
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Figure 4. Same as in Figure[3]but witli (a) 6 = 1.005, and (b) 6 = 1. 
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Figure 5. Same as Figure |3] for = 0.95 (below the critical temperature) but 
calculated with different number of spins and trajectories: (a) N = 10^, Nt = 10^; 
(b) N = 10^, Nt = 10^; (c) N = 10^, Nt = 40 (in fact Nt = 5 suffices). Initial values 



Figures [5] depicts the evolution of the oscillator position toward one of the two 
nonzero equilibrium values for 6 = 0.95, below the critical temperature. To attain 
good agreement between the prediction of the nonlinear equation and the averages over 
trajectories, the number of particles in our simulations has to increase while a small 
number of trajectories (as low as 5) suffices: compare Fig. [5]^a) for = 10^ with|5]^b) for 
= 10^ and with|5]^c) for = 10^. For 6 = 0.95, Figure |5|^c) shows that the nonlinear 
approximation gives again a good description of the oscillator dynamics, accounting for 
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its time evolution to x^q — 0.38, as predicted by (123]) . If the temperature is further 
lowered, the number of spins necessary for the prediction of the nonlinear equation 
to approximate the simulation values again decreases to = lO'^ with Nt = 100 
trajectories. This can be seen in figure [6]^ a) for 9 = 0.9. For 9 = 0.6, figure |6]^b) shows 
that the nonlinear equation predicts a monotonic approach to the equilibrium value 
Xeq — 0.91. On the other hand, the simulation gives underdamped oscillations towards 
Xeq, with a period approximately given by the oscillator natural period. Similar curves 
are found for lower temperatures. The nonlinear equation has the correct equilibrium 
oscillator positions as stable stationary solutions (so it gives the attractors correctly), 
but is not expected to be accurate for temperatures below that given by (15^ . Estimating 
X by its steady value fl28|) . the lowest temperature for which the nonlinear equation flSTl) 
is expected to hold is given by ^ = 0.877. Let us recall that the reason for this is that 
the spin relaxation time diverges as T — )■ [HI [12], [13] , and the separation of time scales 
leading to (I5T]) is no longer valid. 

6. Conclusions 

We have studied a harmonic oscillator subject to a force due to a chain of spins whose 
coupling constant is proportional to the oscillator position. The spins are in contact 
with a thermal bath at constant temperature and evolve following Glauber's dynamics. 
We have shown that the oscillator potential energy is modified by the spins and that it 
experiences a nonlinear friction. The quasi-stationary approximation f[50]) is basically a 
linear theory around equilibrium, which is valid if e <^ 1. Physically, this means that the 
natural oscillator period 2tt/u)q is much larger than the characteristic relaxation time 
of the spins' energy. Then the spins relax to equilibrium over a time scale in which 
the position of the oscillator can be considered roughly constant. The equilibrium 
contribution of Ci accounts for the renormalization of the potential, while the term 
corresponding to the (small) deviation of the Ising system from this "equilibrium" gives 
rise to the friction term. 

The oscillator rest points are the stationary solutions of the corresponding reduced 
dynamical equation. These solutions undergo a supercritical pitchfork bifurcation as the 
bath temperature crosses a critical value. For temperatures above critical, the stable 
equilibrium position of the oscillator is zero, the same as that of the uncoupled oscillator. 
Below the critical temperature, there are two stable symmetric equilibrium positions. 
This pitchfork bifurcation corresponds to a second order phase transition for the 
equilibrium probability of the oscillator-spin system. The oscillator equilibrium position 
is the corresponding order parameter and it plays the same role as the magnetization in 
an effective long range Id Ising system. 

Even when our dynamical equation (IST]) does not give an accurate description of the 
oscillator time evolution for very low temperatures, the equilibrium points are always 
correctly predicted by the solutions of f[25]) . This is not surprising: f[25]) is exact in 
the thermodynamic limit, independently of the value of e, while (15T]) holds only if the 
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Figure 6. Averaged trajectories {x{t)) = x{t) (circles) versus the nonlinear (solid blue 
line) prediction for initial data a;(0) = e^, a;(0) — e and (a) 9 = 0.9 and (b) 9 = 0.6. 
Other parameter values are as in Figure [1] 

characteristic relaxation time of the spins is much smaller than the oscillator natural 
period. This condition is not fulfilled for T — )■ 0, because the spin relaxation time 
diverges in that limit [HI [121 IE] • 
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Appendix A. Normal solution of the system of equations ( I48p - ( 1491) 

The solution of (l48l) satisfying Cq^^ = 1 is 
^(0) ^ ..fx 



Cr =V\ V = tanh ( - ) (A.l) 



7 = TT^ , (A.2) 



2ri 

l + V 

in which Cn^ is bounded as — )■ oo. Inserting these expressions in f H9|) . we get 



^^. + ^.]-2CL^^ = n,-^^, (A.3) 



with the boundary condition C^^^ = 0. This equation can be solved by using standard 
methods for difference equations [21], with the result 

d'^=a„r7", (A.4) 

where 

Then Ci ^ {1 + aie)r] yields dSD]). 
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